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^^ ' Abstract 

Real-time magnetic resonance imaging (MRI) methods generally shorten the measuring time by 

yS^ \ acquiring less data than needed according to the sampling theorem. In order to obtain a proper 

image from such undersampled data, the reconstruction is commonly defined as the solution of an 

S inverse problem, which is regularized by a priori assumptions about the object. While practical 

realizations have hitherto been surprisingly successful, strong assumptions about the continuity of 
image features may affect the temporal fidelity of the estimated images. Here we propose a novel 
^— H ' approach for the reconstruction of serial real-time MRI data which integrates the deformations 

^ , between nearby frames into the data consistency term. The method is not required to be afflne or 

rigid and does not need additional measurements. Moreover, it handles multi-channel MRI data 
by simultaneously determining the image and its coil sensitivity profiles in a nonlinear formulation 
which also adapts to non-Cartesian (e.g., radial) sampling schemes. Experimental results of a motion 
phantom with controlled speed and in vivo measurements of rapid tongue movements demonstrate 
image improvements in preserving temporal fidelity and removing residual artifacts. 
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^ ; 1 Introduction 

Imaging speed is crucial in real-time MRI studies of physiologic processes, ranging from cardiovascu- 
lar imaging to noninvasive monitoring of interventional (surgical) procedures. Because the physical 
acceleration of the data acquisition process is limited by physiologic regulations to prevent peripheral 
nerve stimulation, most strategies ultimately reduce the measuring time by acquiring less data, while 
attempting to preserve the quality of the reconstructed image. A first development along this line is 
the parallel imaging concept which takes advantage of multiple receiver coils to acquire data simultane- 
ously. Such techniques benefit from encoding part of the spatial information of an object into spatially 
complementary coil sensitivities, which are generally unknown and also depend on the actual experi- 
mental condition. Therefore, coil sensitivity profiles are either explicitly pre-calibrated in image space, 
like SENSE [23], or implicitly in /c-space, like SMASH (25] and GRAPPA [9 . Unfortunately, however, 
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such pre-calibration techniques make only suboptimal use of the available data from multiple receiver 
channels, so that respective errors in the estimated coil profiles may lead to artifacts in the iteratively 
optimized image - already for moderate acceleration factors of about two to three. An improved strategy 
is to compute spin density maps and coil profiles at the same time by means of a nonlinear formulation 
of the inverse reconstruction problem [36, 26 . In this case, when aiming for high temporal resolution, 
the use of strongly undersampled data introduces additional ill-posedness to the reconstruction problem. 
In order to stabilize the problem and obtain plausible solutions, it is necessary to incorporate a priori 
information about the unknown object (described by its spin density) and the coil profiles into the 
reconstruction method. In [28] |27J [38], a temporal L 2 -regularization on the object was studied, which 
is, however, usually too weak to remove residual artifacts. For example, temporally flickering artifacts 
are observed for a radial sampling scheme which employs complementary sets of spatial encodings in 
consecutively acquired datasets. In practice, a temporal median filter may effectively diminish residual 
streaking artifacts though at the expense of degrading the true temporal resolution. Alternatively, the 
total variance or total generalized variance were used for regularization [16] , which also reduce streaking 
artifacts, but fail to recover small-scale details of the object and therefore sacrifice spatial resolution. In 
general, regularization methods alone seem to be unable to provide artifact-free images with high spatial 
and temporal resolution in real-time MRI scenarios with pronounced undersampling. 

An alternative strategy for improved image quality is to integrate information about moving object 
features into the reconstruction by exploiting multiple measurements at different (neighboring) time 
points. One of the most effective means for motion compensation in MRI is the 'navigator echo' technique 
and its variants [7] [31] [34] USB ESJ HTJ [T7] [19] . In these methods, a navigator signal is repetitively acquired 
during the scan to extract specific motion information. The need to insert multiple navigator modules 
into the MRI sequence may be avoided by 'self-navigating' techniques which may determine motions from 
the actual data. Early applications [24, 2 used partial or full fc-space data in a block-based or parametric 
manner, but failed to detect complex motions such as elastic deformations. Lately, more flexible motion- 
detection techniques were developed for free-breathing cine MRI studies p~4] [15] [10] [29] , but they can only 
compensate for slow (e.g., respiratory) movements that affect a faster (e.g., cardiac) motion of interest. 
Motion compensation was also combined with conventional parallel MRI reconstructions [20] [33] . 

To overcome the aforementioned limitations, this work presents a novel reconstruction method for 
real-time MRI by integrating the idea of a 'self-navigating' motion into a nonlinear formulation of the 
inverse problem which simultaneously estimates spin density and coil sensitivity profiles. Based on a 
non-parametric motion estimation, the new method generates images with high temporal fidelity and 
reduced residual artifacts. The validation of the approach in comparison to current algorithms employed 
a motion phantom with controlled speed as well as real-time MRI studies of movements during human 
sound production. 

2 Theory 

The proposed method is based on a recently developed reconstruction from highly undersampled radial 
MRI acquisitions with multiple receiver coils [26, 28 . It generalizes the respective data consistency term 
to incorporate an aggregated reconstruction from multiple frames with non-parametric motion correction 
(i.e., AME = Aggregated Motion Estimation) and is schematically outlined in Figure [T] 

2.1 Real-time Magnetic Resonance Imaging 

The real-time MRI data acquisition of the t-th (t G N) frame from multiple receiver coils is given by 

Vt,i = StHpt ■ ct,i) + e t ,i, I G A(:= {1, . . . , TV}), (1) 



where p t denotes the spin density, c t ,i the sensitivity profiles of N individual receiver coils, et,i the 
noise, St : f •->• (f(kt,j))jLi ^ ne sampling operator at the positions (k t j)j / L 1 in fc-space, and T : / i->- 
J f(x)e~ 27rl ^ x,k ^dx the Fourier transform. The goal is to obtain a serial stream of images (pt)teN with 
high spatial and temporal resolution from the measured data (yt,i)teN,ieA- The MRI measuring time per 
frame is mainly determined by the number of samples M (times the physical repetition time TR needed 
for radiofrequency excitation and spatial encoding), which therefore is kept as small as possible. On the 
other hand, for pronounced undersampling conditions, equation (pQ) becomes increasingly ill-conditioned. 
As a consequence, the inversion of the system leads to an amplification of noise which in turn results in 
low-resolution images. Thus, a proper choice of M should be a sensible trade-off between temporal and 
spatial resolution. 

2.2 Aggregated Motion Estimation for Nonlinear Reconstruction 

For explanation purpose, in this subsection, we first assume that the motion (or deformation) (j) t , s from 
p t to pt+s, i.e. (j) t ,s(pt) ~ Pt+s, is known for every s G K, with some KcZ. For example, K can be 
{—2,— 1,0, 1,2}, which corresponds to 5 successive acquisitions. In this case, only two future frames are 
included corresponding to an ignorable waiting time of around 60 ms in our applications. By variable 
substitution, the spin density of the t-th frame satisfies 

St+aFifaAPt) ' c t+s,i) = yt+s,u I e A,s e K. (2) 

Thus, if successive frames rely on complementary data samples in /c-space, the reconstruction takes 
advantage of \K\ • M samples for recovering p t , while the temporal resolution remains unchanged, i.e. 
corresponding to M times the repetition time TR. Accordingly, while keeping the temporal resolution, 
the approach may obtain images with higher spatial resolution from ([2]). For sake of clarity, we rewrite 
([2]) by an abstract nonlinear operator equation 

F t+S (®t : s(x)) = y t+3 , for s e K, with x = ( , / J , Vt+s = (yt+ajheA- 

\ \Ct+s',i)s'eK,ieA J 

Here, F t + 8 : x \-> (St+ 8 F(pt ■ c t+s,l)) leA and $ M : x \-> (0 t , s (/?t), (c t + s ',i)s'eK,ieA) T , for every s e K. 

These equations are solved for the unknown x by a Newton-type method, whose key idea consists 
in repeatedly linearizing the operator equation F t + S (&t,s(%)) = 2/t+s? s £ K, around some approximate 
solution x n , and solving the linearized problems 

Fl +a ($tA x n))$t l8 M( x ~ Xn ) = yt+ s ~ Ft+s(<S>tA x n)), seK, (3) 

As the real-time MRI problem is highly ill-posed, (F/ +s ($^ s (x n ))^>J s (x n )) is not bounded or se- 
riously ill-conditioned. The standard Newton method is not applicable and may not even be well 
defined for noise-free data, because yt+ s ~ F t + S (&t,s(%n)) is not guaranteed to lie in the range of 
Fl+ S ($> t ,s(xn))&t s( x n)- Therefore, some regularization method has to be employed for solving the 
linearized equation (J3j). 

Because only the product of the spin density and coil profiles is determined, the real-time MRI 
problem is undetermined even in the fully sampled case. Although the image may contain fine structures, 
the coil profiles are generally rather smooth. As in [26], this can be ensured by introducing a term 
promoting smoothness which may be given by a Sobolev norm ||/|| m := ||(1 + 6|| • \\ 2 ) m ^ 2 (Ff)(-)\\, 
m G M+ . It penalizes high spatial frequencies by a polynomial of degree m as a function of the distance 
to the centre of /c-space. The object (i.e., its spin density) usually deforms continuously and smoothly 
from frame to frame, so an efficient regularization penalizes the differences between neighboring frames 
to ensure temporal continuity. By combining the standard Newton method and the aforementioned 



regularization, we obtain the well-known iteratively regularized Gauss-Newton method (IRGNM) [TJ [3] 
for solving (|2]) 

h n = argmin{^ \\F( +S (^ s (x n ))^ s (x n )h - (y t+s - F t+s (<f> t:S (x n ))) || 2 



seK 



a n £ || a (i + 6||.|| 2 )^^(c^+/i Ct+s , ! -c° +S) ^' 2 



seK,ieA 

,( n ) _i_ h _ ^0i|2i 






where x n = (p^ n) , (c^ s ^) sG k,^a) T , ft = (h Pt ,(h Ct+sl ) seK) i eA ) T and x = (p£, (c° +s ^) sG k,^a) T the 
initial guess. If the initialization is close enough to the true solution, the choice of a n := ao<? n , with 
< q < 1, is usually sufficient for convergence (cf. [11 J. Because the IRGNM method reduces to 
a Gauss- Newton method for a n = 0, it uses the more robust descent direction at the beginning of 
the iterative process (far from the solution) and the faster convergent algorithm at the end (near the 
solution) . The choice of parameter a serves as a balance between penalization of the spin density and coil 
profiles. The quadratic optimization is equivalent to solve its normal equation, precisely a linear equation 
with a symmetric coefficient matrix. Unfortunately, this linear equation is numerically ill-conditioned 
for large m. A simple preconditioning by the following variable substitution can significantly reduce its 
condition number, making it numerically stable. Let 

A/: I , f* W 



(ct+s,i)seK,ieA J V (Wc t + s ,i)seK,ieA 

with W = J 27 * o a-^l + b\\ • || 2 )- m / 2 . If x := M~ x x, G t+S '•= F t+S o $ M o M, s € K, then an equivalent 
form of IRGNM is given by 

h n = argmin ^ \\G' t+s (x n )h - (y t+s - G t + S (x n ))\\ 2 + a n \\x n + h- x || 2 

h s£K 

^n+1 = %n ~T~ l^n- 

Explicitly, the optimality condition for this quadratic optimization is 



( ^2 G 't+s(xnyG' t+s (x n ) + a n I J h n = ^2 G 't+s(xnT (yt+s - G t+S (x n )) - a n (x n 
\seK J seK 



•So)- (4) 



Its discretized form can efficiently be solved by the conjugate gradient (CG) algorithm. Together 
with the motion estimation described in the next subsection, this strategy represents our novel AME- 
based nonlinear reconstruction method for high temporal and spatial resolution (see Figure [T]). In this 
manner, multiple acquisition frames are exploited for reconstruction with proper motion correction, 
implicitly increasing the sampling rate while preserving temporal sharpness. 

2.3 Motion Estimation 

In general, the differential positional displacements between nearby frames are not known. Therefore, we 
first precompute each frame by the nonlinear inversion (NLINV) method introduced in [26j [28] , which 
corresponds to K = {0} as defined in the previous subsection. Subsequently, the motion is estimated on 
these precomputed images, using the TV-L 1 optical flow model (cf. [32JE]). In detail, the motion (j) t , s 
from p t to pt+s is estimated by (j) t ,s(pt)(%) — Pt( x + u(x)), with u given by the solution to 

min \\p t + Vp t • u - p t+s + v\\i + A|| Vu\\i + /i||Vv||i. 

u.v 
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Figure 1: Flow diagram illustrating the AME-based nonlinear reconstruction method for reconstructing 
the t-th frame, with K = {-2, -1, 0, 1, 2}. 



Because </>t,s(Pt)(x) :== Pt(x -\- u(x)) « p t + Vp*; • u, the auxiliary variable v models the varying 
reconstruction artifacts in different frames. For radial MRI acquisitions, the residual streaking artifacts 
have a relatively low total variation in comparison to the object which contains all local structures. 
Therefore, it is expected that u can only capture the true motion of the object instead of the undesirable 
motion of the artifacts contained in the precomputed images. In order to avoid the impact of outliers on 
motion estimation, we only used the L 1 norm. This non-smooth minimization can efficiently be solved 
by the first-order primal-dual algorithm proposed in [5] . 



2.4 Discretization 

By denoting C t + a ,i(%) := faAP) ' w ^t+a,h with x := (p, (c t+s /^/) s / G K,/'eA) T , for s G If , / G A, a detailed 
formula of (JH is 

E C 't +s ,/(5n)*^^ +s 5 f+s ^C t ' +s ^(5 n ) + a n J | «I„ 



E 

seK,ieA 



C t+s,l(Xn)*F*Sl +s (yt+8,l ~ St+sFCt+s,l{Xn)) - an(Xn ~ Xq). 



This equation will be solved by the CG algorithm which requires repeated application of the operations 
St+ S F and J r *S^ +s . For numerical computation, every function involved needs to be approximated 
by a discretized form of points on a rectangular grid. Since density and coil profiles are compactly 
supported, the Fourier transform T can be computed by fast Fourier transform (FFT) with proper 
periodic extension. 

If the sampling trajectory represents a non-Cartesian radial scheme as used for real-time MRI in 
[26l[28], the computation involving the sampling operators St+ S and 5 t * +s is not straightforward. With 

(St+af,g) = (/, SUsd) and St+af = Ejli ffa+aj), we have S*+ S g = Y<f=i <*(' - k t+s,j)g(h+s,j)- Then, 
J^S^yt+8,1 = Ylf=iyt+s,ie 27ri{kt+s ^^ can be computed by inverse FFT after gridding |22| H2 | H] . or 
nonuniform FFT [8] [T3] . With respect to J r *6'^ +s 5't+ s J r , we have 

(T*S; +s S t+s Tf)(x) = f E <% - kt+sj)e 2iri{x ' v) J f(z) e - 2 ^ k ^^dzdy 



J 



M 



e 



i=i 



>f(z)dz 



= q* f 

where q(x) := ^ 7 =i e 27Vt ^ kt+s ^ x \ It can be computed by two FFTs and one inverse FFT, with Tq given 
by the gridding algorithm. In the Cartesian case Tq equals ones at measured points and zeros elsewhere. 
To sum up, equation (|4]) can numerically be solved in an efficient way. 

3 Methods 

3.1 Data Acquisition 

The proposed reconstruction technique was evaluated for real-time MRI measurements of a motion 
phantom as well as for different parts of the human body in vivo. All studies were conducted on a 3T MRI 
system (Siemens Magnetom TIM Trio, Erlangen, Germany). Continuous data acquisition was achieved 
by using a radiofrequency-spoiled radial FLASH (fast low angle shot) pulse sequence developed for real- 
time MRI [37 . Tl-weighted images were generated by a short repetition time TR (approximately 2 ms) 
and a low flip angle of the RF excitation pulse (5 to 10 degree). A highly undersampled radial /c-space 
encoding scheme was employed with an interleaved arrangement of spokes for five successive datasets 
(i.e., frames). Each single turn corresponded to a full image and contained only a small number of 
spokes that were equally distributed over a full 360°circle in order to homogeneously sample /c-space. To 
prevent aliasing effects from object structures outside the selected field-of-view, a readout oversampling 
by a factor of two was used during data acquisition without compromising imaging speed or signal, for 



Table 1: Acquisition parameters for real-time MRI. 



Acquisition Parameter 


Motion Phantom 


Tongue Movement 


Repetition time (ms) 
Echo time (ms) 


2.28 
1.48 


2.2 

1.4 


Flip angle (°) 
Field-of-view (mm 2 ) 


8 
256 x 256 


5 

192 x 192 


Base resolution 


128 x 128 


128 x 128 


Section thickness (mm) 


5.0 


10.0 


Voxel size (mm 3 ) 
Spokes per frame 


2.0 x 2.0 x 5.0 
9 


1.5 x 1.5 x 10.0 
15 



details see [37]. For human studies, healthy subjects with no known illness were recruited among the 
university students and written informed consent was obtained in all cases prior to each examination. 

The motion phantom consisted of a polyacetal disc rotating with respect to its geometric center. 
Three water-filled tubes with approximately 10 mm diameter were fixed on the disc with a distance 
to the center of 25 mm, 37.75 mm, and 55 mm, respectively. The MRI signals were acquired using a 
32-channel head coil (Siemens Healthcare, Erlangen, Germany) and the measurements were performed 
with three different rotational speeds at angular velocities of 0.5 Hz, 1.0 Hz, and 1.5 Hz, respectively. 

Real-time MRI of the human body was performed in a supine position for studies of the heart and 
movements of the tongue during playing a plastic mouthpiece of a brass instrument. In the latter case 
subjects were asked to perform rapid tongue movements (staccato) at a rate of about 5 Hz. A mid-sagittal 
image was chosen to cover the oropharyngolaryngeal area, while MRI signals were acquired by combining 
a 4-channel small flexible receiver coil (Siemens Healthcare, Erlangen, Germany) and a bilateral 2x4 
array coil (NORAS MRI products, Hoechberg, Germany), using the same setup as previously reported 
for real-time MRI of speech generation [18]. Cardiac MRI was performed during free breathing and 
without synchronization to the electrocardiogram [27] [38] using a 32-channel body coil consisting of 
an anterior and a posterior array with 16 elements each. Online image control employed conventional 
NLINV reconstructions with a post-processing temporal median filter (NLINV-MED). Details of the 
imaging parameters are summarized in Table [1] 

3.2 Image Reconstruction 

All reconstructions were done offline using an in- house software package written in Matlab (R2012a, The 
MathWorks, Natick, MA). In the first step, data from up to 32 receiver channels were combined into 
a small set of 10 virtual channels based on a principal component analysis, as previously described in 
[5] [38]. For the interpolation in /c-space from radial spokes to Cartesian grids, a Kaiser-Bessel window 
function with L = 6, /? = 13.8551 and a 1.5 fold oversampling was used [4 J. To speed up the pro- 
cess, the interpolation coefficients were precalculated and stored in a look-up table. In the next step 
the interpolated data were normalized such that the L 2 norm equaled 100. This allows for choosing 
the reconstruction parameters independent from the data acquisition parameters, which minimizes the 
operator interference and also maintains the quality of the results. 

For AME reconstruction of experimental datasets, we found empirically that it is sufficient to choose 
K = {—2,-1,0,1,2} to exploit the complementary information from 5 successive acquisitions with 
interleaved radial encodings. Numerically, K = {—3,-2,-1,0,1,2,3} gives similar results for both 
simulated and real data (not shown) but increases computational complexity. The preliminary images 
were reconstructed by NLINV, with almost identical parameters for regularization and penalization of 
coil sensitivity profiles. The same initialization was used with the spin density set to ones and coil 
sensitivities to zeros for the first frame. Later both were replaced by the reconstruction results from the 




Figure 2: Principle performance of the motion estimation algorithm (short-axis cardiac MRI). For two 
frames at end diastole (NLINV1) and end systole (NLINV2) motion fields were estimated based on 
their NLINV reconstructions and then applied to NLINV1 as DEFORM. DIFF represents the difference 
between NLINV2 and DEFORM. 

previous frame. For motion estimation, the spatial deformation of the object was implemented with a 
bicubic interpolation, and parameters of the model were set to A = 0.02 and \i — 1.0. 

For comparison, the same data was also reconstructed by the standard NLINV method [271 [38] which 
was implemented in the same software environment. In addition, a temporal median filter was applied 
to the images to reduce residual streaking artifacts as proposed in [28] [27] [38] . It was implemented as a 
post-processing step with a window width covering 5 neighboring images (NLINV-MED). 

4 Results 



4.1 Motion Estimation 

The principle of the proposed motion estimation is demonstrated in Figure [2] using data for the human 
heart. Two frames at end diastole (NLINV1) and end systole (NLINV2) were selected to depict distinct 
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Figure 3: Reconstructions from a phantom (three tubes of water) moving with different constant angular 
speed for (left) NLINV, (middle) NLINV-MED, and (right) AME. The rotating frequencies of the moving 
tubes are (slow) 0.5 Hz, (medium) 1.0 Hz and (fast) 1.5 Hz, respectively. 



differences due to myocardial contraction in preliminary NLINV reconstructions (arrows). The calculated 
deformation of NLINV1 by incorporating the estimated motion is shown as DEFORM. It clearly identifies 
the contraction of the myocardium, whereas the streaking artifacts at the top-left corner of the image 
remain similar as in NLINV1. The example demonstrates that the information of the moving object 
has correctly been captured by the motion estimation, while the image artifacts, which may also change 
from time to time, are appropriately excluded. This can also be visualized in the difference between 
DEFORM and NLINV2 (DIFF), where the artifacts are dominant and the structure of the heart is less 
visible. In the next subsection, we show how the artifacts will be removed, rather than enhanced, by 
aggregating the estimated motions in AME (cf. Figure [1]). 



4.2 AME in Action 

4.2.1 Motion Phantom 

Figure [3] compares reconstructions for a phantom moving at different speeds that were obtained by 
NLINV, NLINV-MED, and the proposed AME method, respectively. For the lowest velocity, all three 
methods produce acceptable results, although the latter two surpass NLINV in reducing streaking ar- 
tifacts (arrow in top row). At moderate velocity, NLINV suffers from stronger artifacts due to faster 
motions, while NLINV-MED even distorts the structure of the fastest moving outermost tube (arrow in 
middle row). The stretched shape of the circular tube is a typical effect from the temporal median filter. 
In contrast, the AME reconstruction offers a proper image with almost no motion or streaking artifacts. 
Finally, for the highest velocity, both NLINV and NLINV-MED result in severely deformed shapes for 
almost all tubes as well as pronounced streaking artifacts. Again, the AME method shows best results 
with only very mild and barely visible artifacts. Furthermore, the signal-to- noise ratio (SNR) of the 
AME reconstruction is higher in all cases compared with the two other methods. 

4.2.2 Human Tongue Movements 

Figure H] demonstrates tongue movements during playing the mouthpiece of a brass instrument. For this 
particular task the tongue tip of the subject had to rapidly move forward and backward touching the 
upper teeth ridge. To better demonstrate the temporal evolution of the motion, a reference line is placed 
at the tongue tip to derive corresponding 2D spatiotemporal intensity profiles. Thus, the flickering of 
the residual artifacts at every 5-th frame is clearly visualized for NLINV. For NLINV-MED the residual 
artifacts are effectively removed at the expense of blurring the tongue movements by the temporal median 
filter. On the contrary, the proposed AME method preserves the sharp intensity changes associated with 
the rapid tongue movements even better than in the original NLINV reconstruction, while at the same 
time successfully minimizing residual streaking artifacts. 

5 Discussion 

In comparison to NLINV reconstructions with and without temporal median filter, the proposed AME 
reconstruction for real-time MRI with pronounced radial undersampling yields serial images with im- 
proved temporal acuity and less residual artifacts. The new approach emerges as an expansion of the 
previously introduced NLINV reconstruction with an aggregated motion estimation which estimates re- 
spective movements from multiple consecutive data sets with complementary spatial encodings. The 
additional information is incorporated into the data consistency term of the nonlinear inverse problem 
for a simultaneous determination of spin density and coil sensitivities. 

Extending other approaches for motion estimation in MRI, the present work is not limited to affine 
or rigid motions. Moreover, the combination of AME with nonlinear reconstruction permits an arbi- 
trary choice of K which defines the set of frames used for reconstructing the actual frame. Extensive 
experimental studies (not shown) demonstrate that a choice of \K\ smaller than the number of frames 
with complementary spatial encodings fails to remove the temporally flickering artifacts. On the other 
hand, choosing \K\ greater than the number of differently encoded frames does not further improve the 
reconstruction but yields comparable image quality. Because the computational complexity increases as 
| If | increases, this behavior explains our choice of K = {—2, —1, 0, 1, 2} in all experiments. As a stopping 
criterion of the iterations in NLINV, NLINV-MED and AME, the well-known Morozov's discrepancy 
principle was initially considered, but it forbids a unique choice of the threshold value for every frame 
because the energy of the signal slightly changes with time even for normalized /c-space data. In our 
applications, we have chosen a fixed number of iterations (i.e., Newton steps) for each method, respec- 
tively, which gives satisfactory results. A data-driven choice certainly appears to be more sensible and 
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Figure 4: Reconstructions of tongue movements in a mid-sagittal plane at a temporal resolution of 33 ms 
for NLINV, NLINV-MED, and the proposed AME. The right column presents spatiotemporal profiles 
(vertical line in upper left image) of the MRI signal intensity for a selected 1.0 s period. 



might be considered in future research. 
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A limitation of the AME method, which may deteriorate its performance, stems from errors in the 
preceding NLINV reconstructions that may lead to an unfaithful motion estimation. This can be seen 
from the slightly blurred temporal profile in Figure |U A natural way to overcome this problem would be 
to run the whole AME procedure at least twice using previous AME (rather than NLINV) reconstructions 
for more accurate motion estimations, though at the expense of further increasing the computational 
demand. In fact, at this time the high computational cost, which is about 20 times that of a comparable 
NLINV implementation on a laptop with MATLAB, is currently the major obstacle for more extended 
practical applications. However, because the computations for each receiver coil and of different frames 
are independent, AME is highly adaptable to parallel computing. Apart from interpolation, the involved 
calculations are simplified to point- wise operations, fast Fourier transform, and scalar products. As is 
shown in Section II-D, the interpolation for non- Cartesian data may be separated from the iterative 
optimization, through a convolution with the point-spread function. These features further ensure a 
possible speed-up by an implementation on graphical processing units. 

6 Conclusion 

This work introduces a new reconstruction method for real-time MRI that offers improved temporal 
fidelity for visualizing rapid dynamic changes. Preliminary results for an experimental phantom and 
in vivo human data demonstrate the practical performance and improved quality which is based on 
the incorporation of estimated object motions into the nonlinear inverse reconstruction process. Future 
improvements are expected by exploiting new regularization methods and by accelerating the computa- 
tional speed. 
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